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ABSTRACT 

We present results from a model of the chemical evolution of protoplanetary disks. In our models 
we directly calculate the changing propagation and penetration of a high energy radiation field with 
Lyman a radiation included. We also explore the effect on our models of including dust grain settling. 
We find that, in agreement with earlier studies, the evolution of dust grains plays a large role in deter- 
mining how deep the UV radiation penetrates into the disk. Significant grain settling at the midplane 
leads to much smaller freeze-out regions and a correspondingly larger molecular layer, which leads to 
an increase in column density for molecular species such as CO, CN and SO. The inclusion of Lyman 
a radiation impacts the disk chemistry through specific species that have large photodissociation cross 
sections at 1216 A. These include HCN, NH3 and CH4, for which the column densities are decreased 
by an order of magnitude or more due to the presence of Lyman a radiation in the UV spectrum. A 
few species, such as CO2 and SO, are enhanced by the presence of Lyman a radiation, but rarely by 
more than a factor of a few. 

Subject headings: astrochemistry, circumstellar matter, ISM: abundances, ISM: molecules, planetary 
systems: protoplanetary disks, stars: pre-main-sequence 



1. INTRODUCTION 

Planets form in protoplanetary disks around low mass 
stars, though the exact mechanism of this formation is 
still not clear. It is therefore important to understand 
the physical conditions in the disks in which they form. 
In particular, studying the chemical history of these disks 
provides us with a wealth of information because the 
chemistry responds to the physical environment. From 
the observed molecular emission we are able to explore 
disk physics such as the overall ionization fraction, which 
impacts on accretion physics, and disk kinematics, a di- 
rect p robe of rotation and turbulence (|Apai fe Lauretta! 
120101 and references therein). In addition, the chemi- 
cal composition of cxtrasolar protoplanetary disks can 
be compared to that of the many remnants from planet 
formation that we observe in our own solar system, such 
as comets, meteorites and asteroids. The composition of 
these bodies provides a chemical memory of the forma- 
tion of our solar system. It is therefore necessary to un- 
derstand the chemical evolution of protoplanetary disks 
so that we can unravel the mystery of how our own solar 
system formed and characterize planet formation else- 
where. 

Because chemical processes respond to the physical 
state of the gas, knowledge of the physical environment 
of the disk is required to properly model the chemistry. 
This includes determining how the chemical composition 
of the disk is affected by the local density and temper- 
ature as well as by radiation from the central star and 
the surrounding environment. Based on observations of 
T Tauri stars, a typical disk system consists of a cen- 
tral star surrounded by a flared dusty disk hundreds of 
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AU in radius. T Tauri stars are low mass young stellar 
objects, masses range from 0.08 - 2 M Q , with emission 
spectra that peak in optical wavelengths, though they 
also e mit in the UV and X-ray regimes (lArdila et al.1 
120021 iKastner et all 11997ft . The optical emission de- 
termines the dust temperat ure within the structure of 
the disk (|Calvet et al.lll991l) . but it is the UV and X- 
ray radiation that provides the ionizing photons needed 
to power molecular chemistry. Due to the reprocess- 
ing of stellar radiation, the disks themselves are flared 
(jKenyon fe Hartmannl ll987). Armed with models of the 
disk physical structure, chemical models have suggested 
that beyond R ~ 40 AU, the chemical structure of the 
disk can be divided into 3 regions: the photon domi- 
nated region, the warm molecular layer and the frozen- 
out midplane (e.g.lAikawa et al.ll200a iWillacv fc Langerl 
l2000t ISemenov et al.l I2004T ). This is in agreement with 
the limits set by observations of molecular emission. Ra- 
diation from the central star as well as the interstellar 
radiation field ionize atoms and dissociate molecules at 
the surface, leading to a photon dominated region at the 
top of the disk where chemical elements exist primarily 
as free atoms. Near the midplane, the temperature is 
low enough that freeze-out of molecules as ices on grains 
can occur. In between is a layer that is warm enough to 
prevent freeze-out, but shielded enough so that molecu- 
lar species can form. The exact locations of these regions 
depend on the chemical species being looked at, due to 
the differences in binding energies and freeze-out tem- 
peratures between different molecules. 

Over the years, chemical models of protoplane- 
tary disks have gro wn increasingly more sophisticated 
( Berg in et al.l [2007, and references therein) . The first 
generation determined the disk composition by assum- 
ing t hermochemical equilibrium in the disk (e.g. iLewisi 
[1971 . Once it was realized that chemical kinetics, 
and not just the solids, played a large role in deter- 
mining the observed abundances, models were devel- 
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oped that direc tly solved the time-dependent reaction 
equat ions (e.g. I Willacy et al.l 119981 : lAikawa fc Herbsd 
1999). The current state-of-the-art models include 
cosmic ray and X-ray ionization, UV radiation from 
the central star, adsorption and desorption onto and 
off of grains, changing dust grain sizes, some simple 
mixin g, surface chemistry and some isotopic chemistry 
(e.g. lAikawa fc Nomural 120081: IWillacv fc Woods! [20091: 
iWoitke et alj|2009t iNomura et al.ll2009l )~ 

One major aspect of protoplanetary disk chemistry 
models is to properly treat the transfer of FUV radiation 
that is incident on the disk. This includes UV radiation, 
both from external sources (typically the Interstellar ra- 
diation field, ISRF) and from the central star (mostly 
the accretion shock, which dominates the stellar radia- 
tion field, SRF, in the UV regime for T Tauri stars) as 
well as X-ray emission from the central star. The two 
different sources of UV radiation have different angles of 
incidence on the disk with the ISRF having a normal an- 
gle of inciden ce, while the SRF has a very shallow angle 
of incidence (| Willacy fc Langer1l2000D . Because of this, 
it has been shown that properly treating the absorption 
and scattering of the SRF is critical in calculating the 
disk chemical evolution as there is little penetration of 
the UV field i nto the disk without scatt ering being taken 
into account ([van Z adclho ff et al"1l2003D . 

In addition to properly treating the UV radiation 
transfer, Lyman a radiation is an i mportant com p onent 
in fu lly de aling with the UV fi eld. iBergin et afl (|2003l 
120041) and iHerczeg et al.l (|2004T l showed that Lyman a 
radiation is dominant in the UV fields of some T Tauri 
stars, and in fact can carry up to 75 % of the UV flux 
from these stars. Additionally, the scattering of the Ly- 
man a radiation in the disk will be different from the 
scattering of the general UV field. Whereas the general 
UV field scatters solely off of the dust grains, the Ly- 
man a radiation scatters from the atomic hydrogen in 
the disk prior to scattering off of dust grains, meaning 
that it needs to be treated sepa rately from the rest of 
the UV field ([Bethel! et alJl201Qf) . The scattering of Ly- 
man a radiation from atomic hydrogen also differs from 
the scattering off of dust grains in that it is isotropic, 
as opposed to scattering off of dust grains which is pref- 
erentially in the forward direction ([Bonilha et al.l [19791 : 
lDrainell2003D . Previous models have included stellar UV 
fields, but no previous model has included an observa- 
tionally motivated radiation field that coupled the reso- 
nant line radiation transfer of Lyman a radiation with 
the full disk chemistry. 

Given the dominance of Lyman a radiation in T Tauri 
stars, the lack of Lyman a radiation in previous disk 
chem istry models is a significant omission (jBergin et al.1 
2003). In particular there are systematic differences in 
photodissociation rates for species that have cross sec- 
tions at 1216 A ([van Disho eck et al.ll2006f) . as compared 
with those that do not, which will dramatically affect 
the chemistry of the disk. For example, ratios such as 
CN/HCN will be affected since CN is dissociated primar- 
ily below 1150 A whi le HCN is sub ject t o dissociation by 
Lyman a radiation ([Bergin et al.l 120031) . The inclusion 
of Lyman a radiation should then increase the CN/HCN 
ratio and p otentially allow for a better fit to the observed 
ratios (e.g. lThi et al][200l . 



Another factor in properly dealing with UV radiation 
propagation is the dust in the disk. Both dust settling 
and dust coagulation will lower the opacities in the upper 
regions of the disk, leading to a flat ter overall disk and 
great er penetration of radiation ([Dullemond fc Nattal 
[2003T ). This increase in penetration depth for the UV 
radiation has been measured observationally by looking 
at the spectral energy distributions of these T Tauri stars 
and measuring the SED slope in the mid-IR wavelengths 
(5-30 fixn) , with more dust s ettling leading to a flatter 
slope (e.g. iFurlan e t al. 2008). These observations show 
that most disks are fit by models with significant dust 
settling included. 

Motivated by these physical effects, we have created 
a protoplanetary disk chemistry model that includes the 
treatmen t of the radia t ion tr ansfer of Lyman a radia- 
tion from lBethell et all (|2010f ) coupled with the full disk 
chemistry to demonstrate how the inclusion of Lyman a 
radiation will affect disk chemistry. We have also looked 
at the effect of dust settling on the chemical composition 
of the disk. In fj2] we describe the model used to cal- 
culate the chemistry in these protoplanetary disks. The 
physical model, as well as the types of chemical reactions 
used, is discussed. In ij3]we present the effects that in- 
cluding dust setting and Lyman a radiation has in terms 
of the chemistry of the disk and in 2] we discuss why 
it is necessary to include them in future disk chemical 
models. 

2. MODEL 

Our numerica l model is based off o f the ALCHEMIC 
code, written bv lSemenov etallpOloh for use on molec- 
ular cloud chemistry, and heavily modified for our pur- 
poses. It is a 1+1D model, which means that while the 
code itself is one-dimensional (ID) over height at a spe- 
cific radius, the code was run at a range of radii to give 
a pseudo-2D result. The disadvantage of such a model 
is that there is no interaction between horizontal zones, 
though with no mixing included this should be a minor 
effect. Our treatment of the UV field, however, was a 
full 2D treatment. The code modeled a flared protoplan- 
etary disk with incident radiation from a central T Tauri 
star as well interstellar cosmic rays. Our model calcu- 
lated the chemical abundances for 639 chemical species 
and 5910 chemical reactions. While most of the reac- 
tions were not time-dependent, we did incorporate a few 
time-dependent reactions into the network which will be 
discussed in more detail below. 

2.1. Physical Model 

In this section we describe the physical structure of the 
disk that was used as well as the radiation field incident 
on the disk. We used the parameters for a typical T 
Tauri star with radius -R* = 2 Rq, mass M* = 0.5 Mq, 
mass accretion ra te M = 10~ s M l r,yr~ 1 and t emperature 
T, = 4000 K (e.g. lKenvon fc Hartmani^flpl. The dust 
composition, taken from lD'Alessio et al.l (|2006f ) . consists 
of two populations, a large grain and a small grain com- 
ponent. In both cases, the grain size distribution is given 

by a power law of the grain radius (a): n(a) = no (jj^) : 
where p = 3.5, a min = 0.005 /im, & max ^ ma u = 0.25 /*m 
and a, max .iarge = f -00 mm. 
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Fig. 1. — Density (left) and Temperature (right) structure of a disk model with e = 0.01 and e = 1. The model was run with the 
parameters for a typical T Tauri star with radius ij» = 2 Rq, mass M, =0.5 Mq, mass accreti on M = 10~ 8 Mpvr~ 1 and temperature T„ 
= 4000 K (e.g. Kcnyon & Hartmann 1995). The dust composition and settling were taken from D'Alcssio ct al. (2006) with a„ 
0.25 /im and a, ma x,large = 1-00 mm. 
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Our model includes the effects of dust settling, param- 
eterized by e, defined as the dust-to-gas mass ratio of 
the small grains in the upper regions of the disk relative 
to th e standard dust-to-gas mass ratio (|D'Alessio et al.l 
2006). A value of e = 1 means that no settling has taken 
place, while a smaller value of e indicates a depletion of 
dust grains in the upper regions of the disk and a corre- 
sponding increase in large dust grains near the midplane. 
We therefore assume that the disk is well shielded in a 
thin layer near the midplane. The general disk temper- 
ature and density structure for dust settling parameters 
of e = 1 and 0.01, shown in Figure [T] are based on the 
assumpt ion that the gas and dust t emperatures are well 
coupled (jKamp fc D ullcmond 2004). When dust settling 
is included, the temperature structure of the disk changes 
which leads to a less flared disk due to the surface density 
being approximated as 



and was 0.02 M for e = 1 and 0.01 M Q for e = 0.1 and 
0.01. 

In order to more clearly understand the gross effects 
of Lyman a radiation on the chemistry of the disk we do 
not include any turbulent mixing in the disk. Addition- 
ally, grain growth in our model is not time dependent, 
but is fixed for each model. While including these time- 
dependent processes in a detailed chemical model might 
be interesting, even with modern computers it is a very 
computationally intensive set of calculations. Many pre- 
vious models of disk chemistry have reproduced results 
in agreement with observations without including these 
processes (e.g.lAikawa fc He rbst 199Et lWillacy fc Langerl 
120001 : lAikawa fc Nomural 120061 ). In comparison, to first 
order diffusion smoothes out abundance variations with 
height and increas e s the depth of the mole c ular layer (e.g. 
Ilgner et al.l 120041: iSemenov et~aT1 120061: iWillacv et all 



20061) . Thus we can still understand the large-scale ef- 
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where Qk is the Ke plerian angular veloci t y, a is a free pa- 
rameter defined by iShakura fc Sunvaevl f 19 730 and c s is 
the sound speed. Because of the dependence of the sound 
speed on temperature, a different temperature structure 
will lead to a different density structure for the disk. 
For the settled models, the temperature of the disk near 
the midplane is increased due to the greater penetration 
depth for radiation in the more settled disks, which will 
play a large role in determining the chemistry of the disk. 
Recent studies have found that protoplanetary disk ob- 
servations are_wejl_^t using dust settling parameters of 
e < 0.01 ([Fur Ian et al.ll2006f ). For the results presented 
in this paper, we compare d e values of 1, 0.1 and 0.01 . 
Based on the work done by I Andrews fc Williams! (|2OO50 
we assume that the inner and outer disk dust popula- 
tions are coupled, and so we treat them identically. The 
mass of the disk depends on the dust settling parameter 



fects of Lyman a radiation on the chemistry without in- 
cluding diffusion or turbulent mixing in our model. 

2.2. UV Field Calculation 

The UV field for the cent ral star that was us ed for 
our model was calculated by iBethell et al.l (|2010l ). who 
explored the radiative transfer and scattering of Lyman 
a radiation in a protoplanetary disk. The UV field was 
calculated such that the Go (100 AU) value, the scaling 
factor between the total flux of the ISRF in the FUV 
band and the total flux in our radiation field at 100 AU 
(|Habind U968). was set to - 700. This value was chosen 
as an average of the Go values observed by iBergin et al.l 
(|2004l ) for a set of T Tauri stars. An example of the UV 
spectra used can be seen in the top panel of Figure [2j 
which shows the UV field for a radius of 208 AU for an 
e = 0.01. The UV field was assumed to be very smooth 
for all wavelengths other than Lyman a, which com- 
pletely dominates the spectra. While there are known 
emission lines in the other portions of the UV spectra 
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Fig. 2. — Top Panel - Sample stellar UV field from lBethell et alj 
( 2010). This example was taken from just below the top of the disk 
for a radius of 208 AU and an e of 0.01. Bottom Panel - The effects 
of scattering and attenuation on the UV field in the disk. Solid lines 
are the continuum, dashed lines are Lyman a radiation. 

(|Bergin et alj |2003). they are orders of magnitude lower 
than the Lyman a radiation and would not significantly 
modify the disk chemistry presented here. The lower 
panel of Figure [5] shows the dependence of the UV field 
on the column density of the disk. Due to geometrical 
dilution, the UV field strength actually peaks slightly in- 
side the disk and not at the surface. The distance from 
the star to a point in the disk is r 2 + z 2 . As z decreases, 
this distance shrinks and the UV field strength increases 
until there is enough material between the star and the 
point for the UV field to be attenuated. 

The Lyman a photon den sity was calculated in the 
two-stage process described in lBethell et al.l (|2010() . Here 
we provide only a high-level description of the aspects of 
the method most relevant to disk chemistry. In order to 
follow continuum FUV and Lyman a photons from the 
star into the disk, the gas and dust opacity must first be 
known. We assume that the pure absorption of both con- 
tinuum and Lyman a photons is simply due to dust. The 
scattering opacity requires more careful consideration 
and is of central importance to the propagation of Lyman 
a. For continuum photons we assume scattering is due 
to dust. For Lyman a the scattering is a combination 
of dust scattering and resonant scattering by H atoms. 
For consistency, we use the same dust grain population 
as that used in iD'Alessio et ail f|2006T ) . The first calcu- 
lation is therefore to establish the distribution of atomic 
H, which is done by balancing the formation of H2 on 
grains with its destruction by UV photons in the Lyman- 
Werner b ands (912-1100AjBlack fc van DishoeckU l987ri. 
Following iSpaans fe Neufeldl (|1997f ) we write the steadv- 
state H2 fraction, /h 2 , as 

"( H 2) = Rn H (2) 
n H ( + 2Rn H ' [ ' 

We adopt the iCazaux fc Tielensl (J2004) expression 
for the grain-surface H2 formation rate R. The sim- 
ple 's hielding function' approach of iDraine fc Bertoldl 
(1996) is used to estimate the H 2 dissociation rate, £, 



ensuring consistency with later chemical computations 
( cj2.3.1[) . The photons are propagated through the (ir- 
regular 2D) spatial grid using st andard Monte Carl o 
radiative transfer techniques (e.g. Be thell et al.1 [2007). 
This approach allows the pote ntially important effects o f 
scattering to be included (e.g. Ivan Zadelhoff et alJ l2003). 
The radiative transfer through the dusty H/H2 distribu- 
tion is iterated until a converged, simultaneous solution 
is achieved for both the H/H2 distribution and UV radia- 
tion field. This iterative approach is typical for problems 
in which the opacity is coupled to the radiation field. It is 
important to note that the H/H2 transition occurs read- 
ily, well above the regions where the bulk of the heavy 
element chemistry takes place. The H/H2 transition de- 
fines the top of the so-called 'warm molecular layer'. In 
this sense the chemistry occurs against a background of 
a fully molecular H2 disk. However, before the Lyman 
a photons can reach these depths they must first pass 
through the upper 'photodissociation layer', composed 
of mostly atomic H. This leads us to the second stage 
of the calculation; the detailed propagation of Lyman a 
photons. We are in essence repeating the UV transport 
calculation, but this time treating the Lyman a propa- 
gation in greater detail. While this two-stage approach 
is not fully self-consistent, it renders the overall problem 
more tractable. 

Although the shape of the stellar Lyman a line bathing 
the disk in T Tauri systems is generally not observed di- 
rectly, it is known to be very wide in the case of TW 
Hya and in other stars with molecule-rich disks (FWHM 
>500kms- 1 . IHerczeg et al|[200l iBeTgin et alj|200l . es- 
pecially when measured in Doppler widths of the disk gas 
(Figure [3]). As a result, the vast majority of Lyman a 
photons scatter off the wings of the resonant-scattering 
Voigt profile presented by HI in the photodissociation 
layer. Even at frequency offsets as large as ±200 Doppler 
widths relative to line center (1216A), the resonant scat- 
tering opacity of H atoms suspended above the disk is 
greater than that due to the accompanying dust (assum- 
ing the typical interstellar dust:gas mass ratio ~ 0.01, 
i.e. e = 1). Of course, this difference is further enhanced 
by the removal of dust via settling. Dust settling also 
enhances the abundance of atomic hydrogen by reduc- 
ing the H2 formation rate (Eqn. [2]). From the point of 
view of a stellar Lyman a photon it is as though the 
disk is covered in an optically thick layer of high-albedo 
material (the atomic hydrogen), which immediately scat- 
ters approximately half of the incident Lyman a photons 
out into space and the remaining photons downward into 
the disk. This downward scattering greatly increases the 
penetrating power of Lyman a photons compared to the 
feebly dust-scattered continuum photons. Due to the 
great width of the Lyman a line, its overall frequency 
evolution as the photons diffuse into the disk is relatively 
unimportant, and the scattering process can be treated 
as conserving photon energy. Nevertheless, in our cal- 
culation we follow the frequency evolution of photons 
using th e detailed theo ry of 'partial frequency redistri- 
bution' (jHummerl [1962). Ultimately these photons are 
lumped together into a bolomctric Lyman a photon den- 
sity, JA(r,z), that is used later in chemical rate equations. 

Because of the proximity of TW Hya, and the lack of 
interstellar gas along the line of site, it is the only system 
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Initial Abundances 3 
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Fig. 3. — The st ellar Lyman a line profile of TW Hya (upper 
panel, taken from IHer czcg ct al. 2002]), and resonant scattering 
(Voigt) profile due to HI atoms in the disk (lower panel). The 
horizontal line in the lower panel represents the dust opacity. 

with detectable Ly a beyond the broad wings. However, 
H2 fluourescent emission is clearly detected in the hand- 
ful of low mass T Tauri disk systems with high resolution 
FUV spect ra. This H2 is pump ed via coincidences with 
Lyman a ( Hcrczcg et a.1.1 l2002f ) . Thus it is likely that 
Lyman a emission is an important facet of the UV field 
in most if not all accreting systems, although the exact 
strength amongst systems needs to be determined. This 
assertion will be quantified with new UV spectra from 
the Cosmic Origins Spectrogr aph on HST, b ut prelimi- 
nary results are in agreement (|Herczegll2010( ). 

In order to determine the effects of Lyman a radiation 
on the chemistry of the disk, we ran two models, one 
with and one without the Lyman a radiation. In the 
case without Lyman a radiation the UV field was iden- 
tical, except that the continuum was interpolated over 
the Lyman a peak to remove it. The results from these 
models will be presented in Sj3l 

2.3. Chemical Model 

For our chemical model we found the abundance of 
each species in the disk by solving rate equations for 
each chemical species at each time-step and zone in the 
disk. The basic equation we solved can be written as 



dn(i) 
dt 



k M n (j) n ( l ) ~ Y k tJ n(j) (3) 



where n(i) is the abundance of species i and fc« is the re- 
action rate for reacting species i and j. At its most basic, 
this equation is a summation of the rates of destruction of 
a species subtracted from the summation of the rates of 
formation. The specific rate equations for each reaction 
used in our network depended on the type of the reaction 
and will be discussed in further detail below. Three-body 



Species 


Abundance 


Species 


Abundance 


H 


1.00 


CO 


1.00E-4 


He 


1.40E-1 


N 2 


1.00E-6 


N 


2.25E-5 


C 


7.00E-7 


CN 


6.00E-8 


NH 3 


8.00E-8 


H 3 + 


1.00E-8 


HCN 


2.00E-8 


S+ 


1.60E-6 


C+ 


1.00E-8 


Si+ 


1.60E-9 


HCO+ 


9.00E-9 


Mg+ 


3.00E-8 


H 2 CO 


8.00E-9 


Fe+ 


2.00E-8 


C 2 H 


8.00E-9 


H 2 0(gr) 


2.50E-4 


OS 


2.00E-9 


GRAIN 


6.00E-12 







a The initial abundances are by number density 
relative to the total Hydrogen abundance 



reactions were not included in our network, as we were 
mostly concerned with the outer regions of the disk and 
three-body reactions become competitive only in the in- 
ner regions of the disk where temperatures and densities 
are very high. Initial abundances of each species were 
take n to be molecular cloud abundances, approximated 
from lAikawa fc HerbstJ (|1999l ) and are listed in Table [U 
Our chemical network was based on the Ohio State 
University Astrophysical C hemistry Group gas-phase 
model from March 200cfl ([Smith et al.l I2004T) . with a 
number of additions that will be elaborated on be- 
low. This network includes some gas-grain interac- 
tion, electron-grain recombination, cosmic ray ioniza- 
tion, cosmic-ray induced photoreactions, ion-molecule re- 
actions, charge exchange reactions, negative ion - neutral 
species reactions, radiative association, associate ejec- 
tion, neutral + neutral — > ion + electron, neutral-neutral 
chemical reactions, neutral-neutral radiative association, 
dissociative recombination, radiative recombination, pos- 
itive ion-negative ion recombination, electron attach- 
ment, photoionization and photodissociation. The re- 
action rates for all of these reactions can be found on the 
OSU website, though only about 10-20% of these reaction 
rates have been measured in the laboratory, which adds 
a certain amount of uncertainty to all chemical network 
calcu lations (e.g. IWakelam et al.l [2006; Vas vunin et al.l 
l200l) . 

2.3.1. Photodissociation 

One addition to the OSU network was the incorpo- 
ration of photodissociation of molecules from measured 
cross sections in the literature. These rates are depen- 
dent on both the strength and the shape of the radiation 
field that is used for the model. Because we were us- 
ing an observationally motivated radiation field, it was 
necessary to properly treat the molecular photodissocia- 
tion. There were two different methods that were used to 
calculate the photodissociation cross sections. For those 
molecules where UV cross sections for the m olecules have 
been measured (|van Dishoeck et al.l 120061 ). the reaction 
rate was determined by 



^photodissociation 



4ttA 



a(X)J x (r,z)d\ s- 1 (4) 



http: / /www. physics. ohio-state.edu/ ~eric/research.html 
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where <r(A) is the wavelength-dependent cross section of 
the molecule and J\(r,z) is the radiatio n fiel d at any 
point in the disk (calculation describe in H2.2\i . For all 
of the species that did not have measured cross sections 
in the literature, the photodissociation rates were calcu- 
lated by 



z (AU) 

280 240 180 120 80 40 20 



^photodissociation 



G Q (100AU) 



ioo.o \ 2 , , ! 



(5) 



(6) 



Go(100AU), as defined above, is the scaling factor be- 
tween the ISRF and the radiation field we used for our 
model and was set to ~ 700. a is the unshielded pho- 
todissociation rate ta ken from the UMIST 2006 database 
fWoodall et al.ll2007i ) . e T accounts for the attenuation of 
the radiation field. Normally, this term would just be 
e~ r , however due to scattering of the radiation the at- 
tenuation needed to be calculated as the radiation field 
for a given height z, (J\, z ), divided by the maximum ra- 
diation field at that radius (J\,z maa: )- 1500 A was chosen 
as a wavelength well beyond that of Lyman a to use for 
this calculation. 

A comparison of the calculated photodissociation rates 
of HCN vs. total column is shown in Figure |U This 
figure illustrates the reaction rate of HCN for three dif- 
ferent methods of calculating the photodissociation rate. 
The dashed line shows the result if we assume a normal 
ISRF scaled to match our radiation field in the FUV, 
Equation (|5]l. but without any scattering included. The 
dotted line includes the e r term which incorporates the 
attenuation and scattering of our UV field at 1500 A 
(jvan Zadelhoff et alJ 12003 Bergin et all [2003). Finally, 
the solid line is the result of using Equation U to calcu- 
late the rate from the wavelength-dependent cross sec- 
tion and our UV field. The spread in the three lines 
at high column densities indicates the importance of in- 
cluding scattering of the UV field in these calculations. 
The difference between the solid and dotted lines occurs 
as a result of using the wavelength dependent cross sec- 
tion instead of just approximating it at 1500 A, which is 
especially important for species with a cross section at 
1216 A that will be affected by the presence of Lyman 
a radiation. Because of geometrical dilution due to the 
flaring of the disk, the maximum UV field is actually a 
few tens of AU below the surface. This effect causes the 
rate calculated from the ISRF alone to be slightly higher 
in the lower column density regions of the disk. 

2.3.2. Self- shielding of CO and H% 

Because the photodissociation rates of CO and H2 are 
dominated by line rather than continuous absorption, 
self-shielding of these two molecules becomes significant 
in the regions where they are optically thick. For our 
model, CO self shielding w as calculated from the rates 
given m iVisser et all (120091) while the H 2 self shielding 
rates were calculated from lLee et ail (j 19961 ). The H 2 pho- 
todissociation rate was 



— This Work 

— FUV approx 

— FUV approx x e T 




Fig. 4. — Plot of three different photodissociation rates of HCN 
at radius R = 208 AU and e = 1. The dashed line (FUV ap- 
prox) assumes a normal ISRF scaled to match the FUV flux of 
a typical T Tauri star without any scattering included. The dot- 
ted line (FUV approx X e T ) includes e r , which incorporates the 
attenuation and scattering of the UV field in the disk. The solid 
line is the result of calculating the photodissociation rate from the 
wavelength-dependent cross section of HCN. 



k H2 ss = 4 x 10- n Go(100yL[/) ( -^£r) 



\R{AU) ) 
x / x e T (r,z,A = 1000A)(s- 1 ) (7) 



with f = 



(life*) for N »i 



> 



and 



f = 1 otherwise. e r is defined in Eqn. [51 though 
in this case it is calculated at 1000 A since that is 
the middle of the range wh ere H 2 is photodissociated 
(|van Dishoeck fc Black) ll988T ). The G (100,4f7) factor is 
used since the equation assumes an ISRF, which needed 

to be scaled up to match the UV fiel d in our mod el. 

For CO, the shielding functions in IVisser et alJ (|2009fl 
were calculated for the unmodified ISRF and thus needed 
to be scaled up to match our UV field. However, there 
was one modification to the scaling technique used to 
calculate the self-shielding rate for H 2 . Because CO does 
not have a cross section coincident with Lyman a ra- 
diation at 1216 A, but the radiation field used for our 
models had a large fraction of its UV flux at 1216 A, us- 
ing a Go value calculated with Lyman a included would 
lead to erroneous results. Instead, the peak at 1216 A 
was removed by linearly interpolating over that portion 
of the curve and that resulting curve was used to calcu- 
late a new Go value used in the CO photodissociation 
calculations. This led to a Go (CO) value of ~150. In 
addition, e T (r,z,A = 960A) was used, with 960 A chosen 
to be in the middle of the CO photodissociation range. 
In both of the self-shielding cases the calculation of the 
rates were time-dependent since they depended on the 
column density of the molecule. 



2.3.3. X-ray ionization and Cosmic-Ray ionization 

Observations of T Tauri stars have shown that 
they emit^ strongly in X-rays, 
erg e 



emit 
tween 10 29 



10 



with luminositi e s be- 
(jFlaccomio et al.l 120091: 
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Feigelson fc Montmerlel 119991: iGlassgold et al.l I1997D 
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Aikawa fc Herbstl (|1999l 120011) discussed the effects that 



this radiation has on protoplanetary disk ionization 
and chemistry, which includes secondary ionization of 
hydrogen gas as well as UV photolysis induced by 
the X-rays. For our model, both of these effects 
have been included u sing the formalisms develope d by 
IGlassgold et~aH (|1997t) and I Aikawa fc Herbstl (|2001f ) and 
using an X-ray luminosity of 10 29 erg s . For the cos- 
mic ray ionization re actions we have follow ed the for- 
malism presented by iSemenov et al.l (|2004f ) with cos- 
mic rays striking the disk vertically on both sides of 
the disk with an attenuatio n column of 96 g cm" 2 
(jUmebavashi fc Na kano 198l]) and an unshielded cosmic- 
ray ionization rate of 1.3 x 10 
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2.3.4. Adsorption & Desorption 

The majority of our chemical network is made up of 
gas species and most of the reactions in our network are 
between two such species. Just looking at gas species, 
however, leaves out a critical component of disk chem- 
istry, reactions between gas species and dust grains. In 
order to self-consistently treat the gas-grain interactions, 
the reaction rates for all such reactions was decreased by 
a factor of e in order to take into account the effect of 
the dust grains settling to the midplane. 

Adsorption is the process by which a gas-phase species 
sticks to the surface of a dust grain and is retained due to 
van der Waals or surface bonding forces. For the low tem- 
peratures that we were looking at in our models we as- 
sumed sticking coefficients of unity ([Burke fc H ollcnbach 
119831 : iBTsschop et al.ll2006D . The reaction rate for species 
adsorption onto grains can then be written as 



8MbT _ x 



(8) 



where a gr is the cross section of the dust grain, approx- 
imated to be ^ r g r - The dust grains were assumed to 
have a radius of r gr = 0.1 fim which leads to a cross sec- 
tion of 3.14 x 10~ 10 cm 2 in Equation ©. This is slightly 
larger than the average dust grain size calculated from 
the distribution in ij2.ll of 0.02 /im due to the fact that 
smaller grains are subject to thermal spikes and it is 
therefore harder fo r molecules to adsorb onto small grains 
(|Leger et al.l ll985). fcs is Boltzmann's constant, /3 is the 
molecular weight of the species, and S is the sticking 
coefficient. 

In our model there were three different methods of des- 
orption, or ways to overcome the binding energy that 
holds a species to the surface of a dust grain after it has 
been adsorbed. 

Thermal Desorption. For thermal desorption, the rate 
of desorption was calculated using the Polyani-Wigner 
relation: 



/ 3.0 x 10 15 k B E b c _ Eh/T 
7r 2 /LtTO# 



(9) 



where Ef, is the binding energy of the molecule and T is 
the dust temperature. Bin ding energies were t aken from 
lHasegawa fc Herbstl (|1993l ) and lWillacvl (|2007l ) with the 



exception of H which wa s set to 1500 K, consistent with 
(|Cazaux fc Tie lcns 2004|). By equating the flux of ther- 
mally desorbing molecules from a grain surface to the 
flux of adsorbing molecules, it is possible to estimate the 
freeze-out temperature for a given species. Doing this, 
iHollenbach et all <|2009f ) estimated freeze-out tempera- 
tures for H 2 of ~100 K, but only -20 K for CO. This 
has two implications for protoplanetary disk chemistry. 
The first is that each species will have a slightly different 
disk structure, in terms of where we see the freeze-out 
region, that is d ependent on the fre eze-out temperature 
of that species (jAikawa et al.l I2002T) . Secondly, the dif- 
ferent temperatures mean that species will freeze-out at 
different locations in the disk, affecting which species are 
available for further chemical reactions in those regions. 

Cosmic-ray Desorption. We included cosmic-ray des- 
orption of ices from dust grains, however this is not a sig- 
nificant desorption mechanism in all regions of the disk. 
As mentioned before, we assumed an attenuation column 
of 96 g cm' 2 for cosmic rays. This means that at radii 
less than ~ 1 AU, where £ is greater than 96 g cm" 2 
even at the midplane, cosmic ray desorption never plays a 
large role in desorption. In contrast, in the outer regions 
of the disk the cosmic ray attenuation is fairly low and 
that, combined with low temperatures at the midplane, 
means that there are regions of the disk where cosmic 
ray desorption is an important desorption mechanism. 
For our model, the cosmic-ra y desorption rate was calcu - 
late d using the formali s m of Hascga wa fc Herbstl ()1993[ ) 
and iBringa fc Johnsonl (j2004f T 

Photodesorption. Photodesorption is the dominant 
form of desorption in the regions of the disk where 
the UV field can reach and the temperature is below 
a species' sublimation temperature. The rate for pho- 
todesorption by UV photons was calculated as: 



k 



'; photodesorption — FjJvY — - Np—N^ 1 (10) 

1* sites ^ice 



where Fjjv is the total UV flux, Y is the yield (the num- 
ber of adsorbed particles ejected per incident photon), 
o-gr is the grain cross section, N S i tes is the number of 
reaction sites on the grain, assumed to be 10 6 , N p =2 
is a correction for the fact that UV photons only pene- 
trate the first few monolayers, n{i) is the abundance of 
a given species, ni ce is the total abundance of all species 
frozen onto grain surfaces and N m is the number of mono- 
layers. This equation is valid for an assumed two or 
fewer monolayers per grain, as laboratory experiments 
by IQberg et al.l (|2007t) indicate that most photodesorp- 
tion occurs in the upper — 2 monolayers. All species 
were assumed to have a yield of 10 -3 , except for CO and 
H2O for which we used recently measured values from 
IQberg et al.l (f2007l . l2009bTi . 

In order to properly treat the Lyman a radiation, the 
normal photodesorption reactions mentioned above were 
calculated with the Lyman a radiation removed from the 
UV field used in Fjjv and instead a separate set of pho- 
todesorption reactions were included for all species that 
had a photodissociation cross section at Lyman a wave- 
lengths. There is little experimental data available on 
photodesorption rates, so we assumed that the physics 
of species on dust grains is similar to the physics of 
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those species in the gas phase. This means that any 
species with a photodissociation cross section at 1216 A 
was assumed to have a related Lyman a photodesorption 
cross section. In order to estimate the photodesorption 
yields of these species, we assumed a constant normal- 
ization factor b etween the gas-pha s e pho todissociation 
cross sections in Ivan Dishoeck et al.l (|2006t ) and the pho- 
todesorption yields. To determine the normalization fac- 
tor we matched the photodissociation cross section for 
H 2 of 1.2 x lCr 17 CTO 2 at A = 1216 A with the pho- 
todesorption yield for H 2 0(gr) of 2.36 x 10~ 3 measured 
bv lOberg et al.l ()200 9b). In symbolic form: 



Yh 2 



x 



OH 2 
OX 



(11) 



where ax is the photodissociation cross section of species 
X and Yx is the photodesorption yield of that same 
species. This new yield was then used in Equation (fTOj) . 
but this time with only the Lyman a radiation flux in- 
cluded. By treating the Lyman a photodesorption sepa- 
rately in this manner, we remove the problem of the UV 
field having most of the flux at Lyman a wavelengths, 
even if the species we are looking at does not have a 
cross section at 1216 A. 

2.3.5. Grain Surface Chemistry 

While the vast majority of our reactions involved gas 
species or molecules adsorbing onto or desorbing off of 
grains, there were two series of grain surface chemistry 
reactions that were included in the model. The first is 
the formation of H2 on grain surfaces. Here we show the 
reaction and the formalism for computing the evolution: 



H + grain H(gr) + grain 



dn 



H{gr) 



dt 

H(gr) + H ^H 2 
if n{H{gr)) > n(gr) : 
dnn 2 
dt 

if n(gr) > n{H{gr)) : 

dnn 2 
dt 



nHn qr < ov > S 



(12) 



n qr nH < ov > S 



n H (gr)nH < ov > S 



(13) 



(14) 



The reaction is calculated in two steps because the lim- 
iting step is for a hydrogen atom to find a grain. Addi- 
tionally, if the abundance of hydrogen atoms on grains, 
H(gr), is greater than the abundance of grains in the 
system (Equation (JT3J) ) , the reaction rate is reduced by 

the ratio J^ 9 ,^ w to take into account the fact that a 

n(H(gr)) 

gas phase atom must collide with a grain in order for 
the reaction to take place. In all cases, the rate coef- 
ficient is < avn > Sh where a is the collisional cross 
section, defined as above as nrg r , vh is the velocity of 
a hydrogen atom assuming a Maxwellian distribution at 
the gas temperature, and Sh is the sticking probability 
of a hydrogen atom, assumed to be one. 
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Fig. 5. — Contour plot of CO and CO(gr) abundances for a model 
with no dust settling included. The dashed lines indicate the r = 
1 surface in the disk. 

The second series of grain surface reactions was oxy- 
gen chemistry t o form water on grains and followed the 
treatment from Hollcn bach et al.r(|2009() . This included: 



O(gr) + H ->• OH(gr) 
H(gr) + O -)• OH(gr) 
H(gr) + OH -> H 2 0(gr) 
OH(gr) + H -> H 2 0(gr) 



(15) 
(16) 
(17) 
(18) 



Laboratory experiments suggest that th ere is rapid for- 
matio n of water ice on grain surfaces (jMivauchi et al.l 
2008), which makes these reactions a necessary addi- 
tion to the chemical network since water is one of the 
dominant oxygen-bearing species. No other reactions on 
grains were included. 

3. RESULTS 

3.1. Basic Model 

For the purposes of the discussion to follow, X(gr) indi- 
cates species X adsorbed onto a grain surface. All abun- 
dances are relative to the total abundance of Hydrogen 
(n(H) + 2 n(H 2 )) by number density. In addition, all of 
the contour plots have been mirrored over the midplane 
to more easily display the features of the disk. Unless 
otherwise stated, all values are taken at a time of 10 6 
years. 

The general structure of our results, demonstrated in 
Figure [5J fits very wel l with the pic t ure d escribed ear- 
lier in ijl] and seen bv lAikawa et al.1 (|2002f ). The fi gurc 
shows + and CO(gr) abundances for a disk with Ly- 
man a radiation and two dust settling values. In the 
upper regions of the disk there is a photodissociation re- 
gion (PDR) where atomic ions and neutrals exist (little 
CO, no CO ice), below this is a molecular layer consist- 
ing mainly of molecules (CO, no CO ice) and near the 
midplane is a cold region where most of the chemistry 
consists of molecules frozen out onto grains (no CO, CO 
ice). This structure can also be seen in the upper-left 
plot of Figure [6l which is a plot of the major carbon- 
bearing species in the disk at a radius of 208 AU vs. 
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Fig. 6. — Plot of the vertical distribution of relative abundances of carbon-bearing (top left), oxygen-bearing (top right), nitrogen-bearing 
(bottom left) and sulfur-bearing (bottom right) molecules at R = 208 AU, e = 1. 



height. In the PDR (z > 200 AU), all of the carbon is in 
the form of C + ions, in the molecular layer (50 AU < z 
< 150 AU) the carbon is in the form of CO, and in the 
freeze-out region (z < 50 AU) the carbon is in the form 
of CO frozen out onto grains. The specific location of 
the CO freeze-out region is sensitive to the temperature 
profile around 20 K, where CO freezes out onto grains. 
As a result, different radiative transfer calculations will 
lead to slightly different loca tions of the freeze-out region 
(e.g. iDullemond eTaI][2002T) . 

Similar plots can be made for the oxygen-bearing, 
nitrogen-bearing and sulfur-bearing chemistry in the disk 
(FigureJS]). For the oxygen-bearing species, the oxygen is 
primarily in the form of neutral O atoms throughout the 
disk until the disk is cold enough for water to adsorb onto 
dust grains. At that point, the oxygen quickly transitions 
from neutral O to H2O ice on grains. There is a small 
layer between these two regimes where molecules such as 
CO and H2O exist, but the molecular layer for oxygen- 
bearing molecules is very narrow. In the nitrogen-bearing 
chemistry, there are again three distinct regions. The 
PDR consists of neutral N which transitions to N2 in 
the molecular layer and then freezes out as N(gr) in the 
freeze-out region. 



The sulfur-bearing chemistry is the most complex. 
While the PDR and freeze-out regions are easy to dis- 
tinguish, consisting of S + and S(gr) respectively, the 
molecular layer is complicated by the presence of mul- 
tiple sulfur-bearing species, such as S02(gr), CS(gr) and 
OCS(gr), that adsorb onto grains at higher temperatures 
than S(gr). 

Focusing on CN and HCN to demonstrate some of the 
more complicated chemical structure, Figure [7] shows the 
abundances of CN and HCN for an e value of 1. Over- 
all we see the general structure of PDR, warm molecular 
layer, freeze-out; however there is also evident substruc- 
ture. The banding feature in these plots occurs due to 
the interaction between multiple formation and destruc- 
tion mechanisms in the model. In the band of lower 
abundances around z = 100 AU, the temperature is low 
enough that there is some freeze-out onto grains, but the 
UV field is still able to penetrate so that there are re- 
actions with PDR species (C+ for HCN and O for CN) 
as well. Lower in the disk, the UV field is unable to 
penetrate as strongly, so the abundance of CN and HCN 
increases briefly until the temperature of the disk is cold 
enough for most of the CN or HCN to adsorb onto grains. 
This is just one example of how the interaction between 
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Fig. 7. — Left Panels - Contour plots of the abundances of CN 
(Top) and HCN (Bottom) for a disk model with e = 1. The dashed 
lines indicate the r = 1 surface in the disk. Right Panel - Plot of 
the vertical distribution of the abundance of CN and HCN in the 
disk at R = 208 AU, e = 1. 

different formation and destruction reactions can com- 
plicate the chemistry of the disk. 

3.2. Dust Settling 

An increase in dust settling will have a large impact 
on th e thermal (Figure [D and chemical s t ructu re of the 
disk (jAikawa fc NomuraM2006l : [Jonkheidl 120061 ) . With 
dust settling included most of the regions of the disk 
near the midplane will have temperatures higher than 
the 20 K needed for CO to adsorb onto dust grains 
( Hasegaw a fc Herbstl 11993ft . Figure [8] shows the abun- 
dance of CO and CO(gr) for e values of 1, 0.1 and 0.01. 
The PDR (top), molecular layer (middle) and freeze-out 
region (midplane) of the disk can be clearly seen in all 
three CO plots. Dust settling increases the size of the 
molecular layer and shifts it towards the midplane at the 
expense of the freeze-out region. This leads to a reduc- 
tion in the overall amount of CO(gr) in the middle of 
the disk. For an e of 1, we see a freeze-out region in the 
bottom third of the disk, while for an e of 0.01, only the 
region at the very midplane of the disk has a significant 
amount of CO(gr). This effect can also be seen in Figure 
[9l which plots C + , the dominant carbon-bearing species 
in the PDR, and CO(gr), the dominant carbon-bearing 
species in the freeze-out region, for e = 1 and 0.01 at a 
radius of 200 AU. Above z > 200 AU the PDRs of the 
two e values are similar, but the molecular layer is much 
larger, and the freeze-out region correspondingly smaller, 
in the case with dust settling included. 

Looking at CO allows us to illustrate the general ef- 
fects of including dust settling in our model. To quantify 
these effects more broadly, we will look at column density 
ratios for a range of representative species at the three 
different dust settling parameters. Figure [TU] compares 
the column densities of species at 250 AU for different 
dust settling parameter values. The majority of species 
plotted, mostly molecular, have larger column densities 
for an e = 0.1 disk than for an unsettled disk with e = 1. 
This is especially true in the outer region of the disk and 
is in agreement with the results in Figure [5] since the 
more settled disk has a larger molecular layer. 

Analyzing the results in more detail at R = 250 AU, 
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Fig. 8. — Contour plot of of CO and CO(gr) abundances for 
models with 3 different e values. The dashed lines indicate the r 
= 1 surface in the disk. 
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Fig. 9.— Plot of C+ and CO(gr) vs. height at R = 208 AU for 
dust settling parameters of e = 1 and e = 0.01. 



the model with the most extreme dust settling, e = 0.01, 
shows more variations from the unsettled disk than the 
e = 0.1 model does. This is unsurprising, as more dust 
settling will allow the radiation to penetrate further into 
the disk and affect the chemical composition more. Look- 
ing at the nitrogen-bearing species specifically, HC3N, 
N2H" 1 ", NH and NH3 are all significantly depleted in the 
settled disk as compared with the disk with no dust set- 
tling. In comparison, CN is slightly enhanced. Almost 
all of these differences can be traced back to the increase 
in photodissociation rates in the settled disk. HC3N and 
NH 3 are directly photodissociated by the UV radiation, 
while NH is destroyed by reactions with neutral O that 
is produced by CO and O2 photodissociation. In com- 
parison, N2H+ is destroyed through reactions with CO, 
which is enhanced due to the wider molecular layer in 
the settled disks. Because the settled disk is warmer 
near the midplane, most of the nitrogen-bearing species 
are unable to freeze-out onto dust grains, which leads to 
an excess reservoir of N in the gas phase. This also serves 
to increase the column density of CN, which is formed 



11 



N & S Species 



Column Density Ratios - 250 AU 



C & O Species 



HCN 



# £■ <? _d> * . 



■ c j ^ ^ ^. ^ c? o -. g * p 



Fig. 10. — Column density ratios for models with different dust 
settling parameters, computed at 250 AU. 



by N reacting with CH and C2- Surprisingly, HCN does 
not show a significant depletion in the highly settled disk, 
despite being strongly dissociated by UV radiation. This 
is due to the formation of HCN from neutral N working 
to counteract the photodissociation. 

Sulfur-bearing species show a similar pattern, with the 
higher temperatures in the settled disk leading to very 
few species frozen-out onto grains. The result of this is a 
significant increase in the column density of species such 
as S, S + and SO in the gas-phase. SO2 is strongly pho- 
todissociated, adding to the enhancement seen in SO and 
O. For the carbon-bearing species, C and C + are both 
slightly enhanced due to the lack of frozen-out species 
as well as photodissociation (of CO and C respectively). 
C2H4 and CH4 are both photodissociated by the UV field 
as expected, leading to the significant depletion of these 
species in the settled disks. 

The lack of water frozen out onto grains in the settled 
disks due to photodesorption leads to an increase in the 
column density of gas-phase oxygen-bearing species such 
as O. The column density of OH is not enhanced as it 
it destroyed due to reactions with S + , which is also en- 
hanced in the settled disk, while H2O is destroyed due 
to photodissociation. For some species, such as CO2 and 
H2CO, an increase in the photodissociation leads to a 
d epletion in the column den sity. 

lAikawa fc Nomural ([2006) explored the vertical distri- 
butions of molecules in disks with larger grains and found 
that the peak abundances are reached at higher A„ val- 
ues for the cases with larger grains. Our results, (e.g. 
Figure 1101) are in general agreement with this finding. 
They found that the column densities of most molecular 
species did not depend on the size of the dust grains, with 
the exception of HCO + and H3" 1 ". The HCO + abundance 
was lower at higher densities, such as disks with larger 
grains and H3 + was found to be abundant in the cold 
midplane, which was thinner in the models with dust 
settling. Our models found a similar trend for H3 + , with 
the abundance decreasing as dust settling increased and 
the cold midplane shrank. HCO + also showed a decrease 
in column density as dust settling was increased. 
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Fig. 11. — Top Panels - Contour plot of the abundance of CN 
(Left) and HCN (Right) for a model with e = 0.01 and no Lyman 
a radiation. Bottom Panels - Contour plot of the abundance of 
CN (Left) and HCN (Right) for a model with e = 0.01 and Lyman 
a radiation included. 

3.3. Lyman a 

Largely depending on whether they have a photodisso- 
ciation cross section at 1216 A or not, some species are 
dramatically affected by the presence of Lyman a radia- 
tion while others are largely unaffected. CO falls in the 
latter category, with no cross section at 1216 A and very 
little change between the models run with and without 
Lyman a radiation included. A clear example of the in- 
teraction between Lyman a radiation and a cross section 
at 1216 A can be seen in Figure QT] CN does not have 
a photodissociation cross section at 1216 A and, as a re- 
sult, we see very little change in abundance between the 
models with and without Lyman a radiation included. 
HCN, in comparison, has a large cross section at 1216 
A and is strongly depleted when Lyman a radiation is 
included in the model. 

However, because of the numerous interconnected reac- 
tions, solely knowing the photodissociation cross section 
is not enough to determine the effect that adding Ly- 
man a radiation will have on the abundance of a species. 
H 2 and OH both have photodissociation cross sections 
at 1216 A, but despite this neither species shows a signifi- 
cant depletion due to the presence of Lyman a radiation. 
Figure [12] shows that OH is actually slightly enhanced 
and H2O is slightly depleted when Lyman a radiation is 
included in the model. These responses are because of 
the importance of Lyman a photodesorption. H20(gr) 
will photodesorb into both H2O and OH, mitigating the 
depletion of these species that would otherwise be seen 
due to photodissociation. Because H20(gr) is the dom- 
inant oxygen-bearing species for the regions where the 
temperature is below 100 K, this photodesorption plays 
a large role in determining the predicted abundances of 
OH and H2O. Our treatment of photodesorption seems 
reasonable, as our maximum abundance for H2O of a few 
ti mes 10~ 6 is only sl i ghtly higher than the results found 
in iHollenbach et all (|2009ft and iDominik et~ail (|2005lh 
Assuming that the H 2 abundance is determined solely 
by a balance between photodesorption and photodisso- 
ciation, a reasonable assumption based on the rates of 
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Fig. 12. — Top Panels - Contour plot of the abundance of OH 
(Left) and H2O (Right) for a model with e = 0.01 and no Lyman 
a radiation. Bottom Panels - Contour plot of the abundance of 
OH (Left) and H2O (Right) for a model with e = 0.01 and Lyman 
a radiation included. 

those reactions compared to the other formation and de- 
struction reactions of H2O, we can calculate an estimated 
abundance of H2 O in the disk by setting those two reac- 
tion rates equal: 



n(H 2 0) x kphotodissodation (Eqn. g}) = 

n(H 2 0(gr) x k ph otodesorption {Eqn. <jT0j)) 



n(H 2 0) / a{X)J x d\ 



F UV Y- 



N. 



^ f «^0 W ) (19) 



Here Fjjv = J J\d\. Using an average cross section (a) 
for H2O of 5 x 10 _18 cm~ 2 , a photodesorption yield (Y) 
of 2.36 x 10 _3 and assuming that the abundance of all 
non H2O species frozen onto grains is negligible, rii ce ~ 
n H 2 0(gr)i we calculate a maximum H2O abundance of 
~ 10~ 6 for the case of e = 0.1, which matches our results. 

To more generally see the impact of Lyman a radiation 
on the chemistry of the disk, we have plotted the ratio of 
column densities for models with and without Lyman a 
radiation included for a variety of species in Figures 1131 
In general, the Lyman a radiation has a more dramatic 
effect when dust settling is included in the model. The 
more settled disks allow greater penetration of the Ly- 
man a radiation, thereby affecting chemistry in a much 
larger region of the disk. In the case with no dust set- 
tling, since the gas and dust opacities are connected the 
Lyman a photons are shielded before they can have much 
of an effect on the chemical abundances. Unless other- 
wise stated, the rest of the analysis in this section will be 
discussing the e = 0.01 the larger effects there 

are easiest to see and because recent observations indi- 
cate that many protoplanet ary disks have significant dust 
settling (|Furlan et al.ll2008D . 

Looking at the carbon-bearing species, only a couple 
of species are strongly depleted. Species, such as C2H4 
and CH4, are significantly dissociated in the presence of 
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Fig. 13. — Column density ratios for models with and without 
Lyman a radiation included, computed at 250 AU. 



Lyman a radiation. Other species that we would expect 
to be destroyed, such as C2H2, are not. This lack of 
depletion is mostly due to an increase in the gas-phase 
abundance of these species that is caused by Lyman a 
photodesorption, similar to what happens with H2O. For 
these species the photodesorption and photodissociation 
mostly balance each other out, leading to only small 
changes in the column densities of those molecules. CO2 
is slightly enhanced due to reactions with O and OH, 
both of which are enhanced in the disk, which counter- 
acts the CO2 photodissociation. 

The nitrogen-bearing chemistry is similar to the 
carbon-bearing species. A couple of species are photodis- 
sociated by the Lyman a radiation as expected, HC3N, 
HCN and NH3, but most of the chemistry is largely un- 
affected. 

The sulfur-bearing species show the most dramatic dif- 
ferences in column densities when Lyman a radiation is 
included. SO2, as expected, is significantly depleted in 
the presence of Lyman a radiation due to its large pho- 
todissociation rate, though only in the highly settled disk 
with e = 0.01. In the less settled disks, formation re- 
actions between SO and OH or O are quick enough to 
counteract the photodissociation. SO also has a large 
cross section at 1216 A, but is actually enhanced by the 
presence of Lyman a radiation. This occurs largely due 
to Lyman a photodesorption of SO(gr), the photodisso- 
ciation of SO2 and the numerous reactions with O and S 
atoms, both of which are enhanced in the presence of Ly- 
man a radiation. For OCS and S, the slight enhancement 
in abundances are due to Lyman a photodesorption of 
SO which is then photodissociated to form other sulfur- 
bearing species. 

Due to the fact that we have very few on-grain sur- 
face reactions in our network, we expect to see a correla- 
tion between the gas-phase and on-grain abundances of 
species. Specifically, species that are photodissociated in 
the gas-phase will be destroyed before they have a chance 
to freeze out onto grains. Additionally, due to the way 
we dealt with photodesorption, even if those species that 
are sensitive to Lyman a radiation did adsorb onto dust 
grains they would be quickly photodesorbed. As a re- 
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Fig. 14. — Column density ratios of on-grain species for models 
with and without Lyman a radiation included, computed at 250 
AU. 

suit, the on-grain species that are strongly affected by the 
presence of Lyman a radiation are those that have gas 
phase photodissociation cross sections at 1216 A, with 
the effects being even more dramatic than for the equiv- 
alent gas-phase species. As can be seen in Figure [TH 
which plots the column density ratios of some important 
on-grain species at 250 AU, species such as C2H4, CH4, 
HCN, NH 3 , SO and S0 2 are all strongly depleted in the 
model with e = 0.01, where Lyman a radiation is able 
to penetrate to the freeze-out region. These species are 
either destroyed in the gas phase before they can freeze- 
out on the grains or they are photodesorbed quickly once 
on a grain surface. One caveat to these results is the 
potential for hydrogenation of molecules on grains. In 
our models we have accounted for the hydrogenation of 
oxygen, which is the main mechanism for water ice for- 
mation. This water can be photodesorbed at the disk 
surface producing water vapor in the molecular layer. It 
is possible that grain surface hydrogenation of nitrogen 
or carbon could perhaps lead to an increased formation 
of species such as NH 3 (gr) and CH 4 (gr), though this is 
not included in the current model. 

Our initial expectations were that the abundance 
changes that occurred would be strongly correlated with 
the size of a species' photodissociation cross section at 
1216 A. This is clearly not the case and in fact there 
seems to be little to no correlation at all between the two. 
The reason for this is that the interdependence of species 
matters more than just the cross section of a species at 
1216 A. For example, the carbon chemistry in the disk 
is completely dominated by the reactions involving CO, 
as seen in Figure [51 In the PDR at the top of the disk, 
all of the carbon is in the form of C + . Further down 
in the disk where there is not enough radiation to pho- 
todissociate the CO molecule, the carbon exists in the 
form of CO. Finally, at the midplane the temperatures 
are low enough for CO to freeze out and the carbon is in 
the form of CO ice on the grains. There is an interesting 
region just between the molecular and freeze-out regions 
in which the temperature is too high for CO to freeze 
out, but cold enough for other molecules to adsorb onto 



dust grains. In this region we start to see the presence of 
more complex molecules frozen out onto the grains, such 
as C3, CH4, HCN and C2H2. Even though the reactions 
that form these molecules are fairly slow compared with 
the rest of the carbon chemistry, since they are being 
depleted from the gas phase through adsorption there is 
enough time to build up a substantial abundance of these 
more complex species and remove that carbon from the 
rest of the network. 

CO does not have a photodissociation cross section at 
1216 A, which means that any species connected to CO 
will have only minor changes in abundances when Lyman 
a radiation is included. This effect can be seen clearly at 
250 AU for species such as HCO + , which is closely con- 
nected to CO, and CH4, which is more weakly connected. 
The HCO + abundance is largely unaffected by the pres- 
ence of Lyman a radiation, while the CH4 is strongly 
depleted. 

The oxygen-bearing chemistry in the disk is actually 
fairly simple in that, like the carbon-bearing chemistry, 
it is dominated by a single species, H2O. As seen in Fig- 
ure [5J the PDR at the top of the disk is made up com- 
pletely of free O atoms, which freeze out as water ice as 
soon as the disk is cold enough to allow it, around 100 K. 
There is a little oxygen in CO and CO(gr) in the molecu- 
lar and freeze-out regions respectively, but these are less 
abundant by a factor of a few than the O and H20(gr). 

The nitrogen-bearing chemistry is similarly straight- 
forward. All of the nitrogen in our model starts in the 
form of N2, but is quickly dissociated to form neutral N 
in the PDR. Neither N2 nor N has a large photodissoci- 
ation cross section at 1216 A, so this process is largely 
unaffected by the presence of Lyman a radiation in the 
UV field. In the molecular layer, a range of nitrogen- 
bearing species exist, including NH3 and HCN. Both of 
these have strong photodissociation cross sections at 1216 
A, so they are destroyed in the highly settled disk models 
leading instead to a molecular layer made up primarily of 
N2. Near the midplane, the nitrogen exists in the form 
of N(gr) when the temperature is low enough for it to 
freeze out, which can be seen in Figure [51 

The sulfur chemistry is a little more complicated. Fig- 
ure [6] indicates that the primary sulfur-bearing species 
are S + in the PDR and S(gr) in the midplane. However, 
the sulfur species with large cross-sections at 1216 A are 
SO and SO2. As was discussed above, the photodissoci- 
ation of SO2 leads to the formation of SO, which is also 
formed through photodesorption and reactions between 
O and S atoms. This means that despite the large cross 
section at 1216 A, SO is actually enhanced in the pres- 
ence of Lyman a radiation, though the SO(gr) abundance 
is depleted. 

3.4. Comparison with Observations 

We have primarily focused on the outer parts of the 
protoplanetary disk in this paper because that is the re- 
gion of the disk probed by current millimeter and sub- 
millimeter observations. In Table [2] we compare the re- 
sults of our model with observations taken of three disks: 
DM Tau, LkCa 15 and TW Hya. The column densities 
listed for our model were taken at a radius of 250 AU and 
a time of 10 6 years for each of the dust settling parame- 
ters. For each species in the table, the top line includes 
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Fig. 15. — Figures comparing our model with e = 0.01 (except for DM Tau, which is e = 0.1) and observations of three disks . The gray 
bars and open circles represent the observations, with error bars of 50%. The triangles show our model with Lyman a radiation included 
and the squares sh ow our model with Lyman a radiation not included. 
top left - DM Tau HAikawa et alj|2003 ), t = 0.1. 

top right - LkCa 15, Interferometric observations (Qi 2001; O.i ct al. 2003; Aikawa ct al. 2003). 
bottom left - LkCa 15, singl e-dish observations jT hi ct al. 2 0041). 
bottom right - TW Hya <Thi et al.ll200l Ivan Dishoeck et alj|2003ft . 



Lyman a radiation while the second line does not. In 
addition, we present figures of the best fit models vs. 
the observations with and without Lyman a radiation. 
These can be seen in Figure [T5l 

We did not produce source-specific models for each 
of these objects, but instead show our results from the 
generic T Tauri star model described in §2.11 The ob- 
jects listed here are all transition disks, which gener- 
ally have lower ac cretion rates and larger disk masses 
(|Naiita et al.l l2007). Despite this, we feel that our mod- 
els can still provide a decent approximation to the obser- 
vations for two reasons. First, M/a, which determines 
the density structure, is roughly consistent between our 
model and these objects. For our model, M = l.Ox 10~ 8 
anda= l.Ox 1CT 2 so M/a = l.OxlCT 6 . Second, the col- 
umn densities derived from observations are very sensi- 
tive to assumptions about the physical structure, adding 
a level of uncertainty to those values. For most of the 



species shown here the addition of Lyman a radiation 
does not affect the column densities much. However, for 
the species that are affected by Lyman a radiation, HCN 
and HNC, the models with Lyman a radiation included 
often produce a better agreement with the observations 
than those without it. The best fit model to each set of 
observations was determined using a two-sample ks-test 
of the observations with each of the 6 models (3 dust val- 
ues, with and without Lyman a radiation). In all cases 
the statistics are severely limited as the best fit is deter- 
mined by so f ew data point s. 

DM Tau - iCalvet et all (|2005[) determined that this 
object has a stellar mass of 0.65 M Q , a stellar radius of 
1.2 Rq, a stellar temperature of 3720 K, an accretion 
rate of 2 x 10~ 9 M yr^ 1 and an a of 0.0015, leading to 
an M/a of 1.33 x 10~ 6 . It has a disk mass of 0.05 M Q 
and the S EP was best fit with e = 0.1. The observations 
are from iDutrev et alj (|1997f) who determined the gas 
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TABLE 2 

Comparison of the Model with Observations 



o pec ICS 


6 


- i 


Models 
e = 0.1 
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01 
U.U-L 


LJ 1V1 IdU 


LkCa 15 
Interferometer 15 Single dish c 


± vv Jnyd 


H 2 


8.3 


(22) 


3.8 


(22) 


3.1 


(22) 


3.8 (21) 










8.3 


(22) 


3.8 


(22) 


3.1 


(22) 










CO 


1.3 


(18) 


6.2 


(18) 


5.1 


(18) 


5.7 (16) 


1.7 (18) d 


1.9 (16) 


3.2 (16) 




1.3 


(18) 
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(18) 


5.1 
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(11) 


2.1 (12) 


2.4 (13) 


1.8 (12) 


8.5 (12) c 




2.5 
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7.4 


(13) 
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(13) 
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1.4 


(12) 


7.3 


(12) 


8.3 


(12) 


9.1 (11) 


<5.4 (12) 




<1.4 (12) 




2.6 


(12) 


4.6 


(13) 


4.7 


(13) 




CN 


5.6 


(12) 


1.5 


(13) 
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(12) 


9.5-12 (12) 


9.7-25 (13) 


1.5 (13) 


6.6 (13) 
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(13) 


9.8 
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H 2 CO 


9.6 
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1.3 


(13) 


1.1 


(11) 


7.6-19 (11) 


7.1-51 (11) 


<8.0 (11) 
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(11) 


1.0 


(13) 
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(11) 


HCO+ 


3.2 


(12) 


6.1 


(11) 


3.2 


(10) 


4.6-28 (11) 


1.5 (13) 


3.3 (11) 


8.5 (12) c 




3.2 


(12) 


7.0 


(11) 
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C 2 H 
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(12) 
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(13) 


1.3 


(13) 


4.2 (13) 










6.9 


(12) 


2.4 


(13) 


1.3 


(13) 










N 2 H+ 


3.5 


(11) 


3.1 


(11) 


1.1 


(10) 


<7.6 (11) 


3.1 (13) d 


<1.4 (12) 


<1.0 (13) 




3.5 


(11) 


2.6 


(11) 


1.0 


(10) 











Note. — Model column densities (cm -2 ) are calculated at t = 1 x 10 6 yr, R = 250 AU. For 
each species, the first row inclu des Lyman a ra diation while the second line does not. The table 
of observations was compiled by I Willacvl {2007). 

HAikawa eFaLl i|2002l . derived from lDutrev et al.lfl997l ) 
b IQil 1120011) unless otherwise noted 
c IThi e t al. (2004]) unless otherwise noted 

d lQi et all (12001) 

e IvanD ishocck etaTI IpOOoj ) 
f IAikawa et all l|2003l ) 



density distribution of the disk and used this to derive 
the average fractional abundances with respect to H2, 
assuming that the fracti onal abundanc e s wer e the same 
everywhere in the disk. lAikawa et al. (2002) used this 
data and integrated vertically to get the observed col- 
umn densities listed in Table [2j Our models for e = 0.1 
produced the best fit to the data, with the model with- 
out Lyman a radiation included actually fitting the data 
slightly better, though the ks-test for the models with 
and without Lyman a radiation were similar (Figure 1151 
top left). 

LkCa 15 - lEspaillat et all (|2007l ) found that this ob- 
ject has a stellar mass of 1.1 M , a stellar radius of 1.7 
R,0, a stellar temperature of 4350 K, an accretion rate 
of 2.4 x 10" 9 M yr" 1 and an a of 0.0006, leading to 

an M/a of 4.0 x 10~ 6 . It has a disk mass of 0.1 M 
and the best fit SED model was with e = 0.001. Look- 
ing at our results for e = 0.01 we find that our models 
are actually a slightly poorer fit to the interferometric 
data when Lyman a is included (Figure [T5l top right), 
but our models with Lyman a radiation included and 
e = 0.01 produced the best fit to the single dish observa- 
tions (Figure IT5| bottom left). One of the largest differ- 
ences between our model and the single-dish data occurs 
in the column densities of H2CO, where our models sig- 
nificantly underestimate the observed values. This is not 
unexpected, as lOberg et al.l (|2009af ) has suggested that 
on-grain reactions can be a critical formation pathway 
for H2CO, which is not currently included in our model. 

TW Hya - This object has been studied extensively 
and has been found to have a stellar mass of 0.6 M , 
a stellar radius of 1.0 Rq, a stellar temperature of 4000 



K (jCalvet et all 120021), a n accretion rate of 3 5 xlO" 9 
Mq yr" 1 (|Espamatll2009l ) and an a of 0.0025 (jEspaillatl 



120091) . leading to a n M/a of 1.4 x 10 ~ 6 . It has a disk 
mass of 0.06 M (|Calvet et all I2002D and th e best fit 
SED model was with e = 0.01 (jEspaillatl I2009D . For the 
observed column densities, the disk was assumed to be 
165 AU in radius. Our model with e = 0.01 and no 
Lyman a radiation included produced the best fit to the 
observations, though again the fit for the model with and 
without Lyman a radiation were very similar (Figure [151 
bottom right). In this case the inclusion of Lyman a 
radiation actually produces too little HCN in the model. 

Column densities for a larger subset of species from our 
model can be seen in Table G2 for e = 0.01, R = 250 AU 
and time steps of 3 x 10 5 and 1 x 10 6 years. 

4. DISCUSSION & SUMMARY 

The existing datasets are quite sparse and our ability 
to compare between model and observations is somewhat 
hampered by the methods chosen for column density de- 
termination. With these caveats in mind, within our 
models the best match is obtained with both dust set- 
tling and Lyman a radiation included. Since both of 
these features are known to be present in accreting gas- 
rich disk systems, in particular TW Hya, DM Tau, GM 
Aur, and LkCal5, we believe that these physical effects 
need to be incorporated into protoplanetary disk chem- 
istry models. The dust settling had the expected effect 
of incre asing the UV penetration as t he dust settling in- 
creases (|Dullemond fc Domi nik 2004). This changed the 
size and location of the molecular layer and freeze-out 
regions, which in turn affected the column densities of 
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TABLE 3 
Column Densities at e = 0.01, R = 
250 AU a 



Species 3 X 10 5 yrs 1 X 10 6 yrs 
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a Column Densities are all in cm 2 

molecules in the disk. Most of the differences that oc- 
curred between the different dust settling models came 
about due to an increase in photodissociation reaction 
rates as dust settling increased. In addition, the smaller 
freeze-out regions of the settled disks allowed species such 
as N 2 H + to be depleted due to the increased molecular 
layers in those disks. Some of these changes were multiple 
orders of magnitude effects and since most disks appear 
to be heavily settled (|D'Alessio et al.ll2006t iFurlan et all 
120061) leaving dust settling out of protoplanetary disk 
models will produce results that are less likely to match 
observations. For the few objects where we were able 
to compare our results with observations, all three disks 
were fit best with a disk where some dust settling was 
included. 

As with dust settling, the inclusion of Lyman a radi- 
ation transfer to the protoplanetary disk model had a 
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significant effect on the resulting chemical abundances. 
While not all species with large photodissociation cross 
sections at 1216 A were heavily depleted, many species 
had column densities that were changed by an order of 
magnitude or more. Some of the most depleted species 
included HC 3 N, HCN, NH 3 , S0 2 , C 2 H 4 and CH 4 , all of 
which were strongly photodissociated by the inclusion of 
Lyman a radiation. In addition, some species had col- 
umn densities that were enhanced, such as S, SO and 
CO2, due to the photodissociation of other species with 
cross sections at 1216 A. 

For some species, the effects of Lyman a induced pho- 
todesorption played a large role in determining how that 
species was affected by the inclusion of Lyman a radia- 
tion. Species such as H2O and OH had their abundances 
increased due to Lyman a photodesorption which coun- 
teracted the destruction of these molecules by photodis- 
sociation. 

Comparing the results from our models with the ob- 
servations of three disk systems indicated that in general 
the models with Lyman a radiation included produced 
similar or slightly better fits to the data. However the 
statistics on these fits were very limited as only a hand- 
ful of species were observed and only two of those were 
species that were affected by the presence of Lyman a 
radiation. Despite this, we feel it is promising that the 
models with Lyman a radiation included in them often 
produced better fits than those without Lyman a radia- 
tion. 

As mentioned previously, the most molecule rich disks 
all have indications of either strong Lyman a emission 
(e.g. TW Hya) or broad Lyman a wing emission (e.g. 
DM Tau, GM Aur, LkCa 15). In both of these cases there 
is also indirect evidence for Lyman a emission through 
fluorescent lines of molecular hydrogen. Therefore we 
believe that Lyman a emission is an intrinsic feature 
of most, if not all, accreting disks. Our models have 
shown that the inclusion of both Lyman a radiation (and 
specifically Lyman a radiation transfer) and dust settling 
are important components for understanding disk chem- 
istry. At present the observational data do not allow 
for definitive statements regarding which physical effect 
dominates, but this will change as we head into the age 
of ALMA. 
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